Gene expression of rabinbow mice

Differential expression

load("DataSummary/DESeq_Cohort1.RData")
options(stringsAsFactors = FALSE)

mTable = DESeq_Cohort1$rList[[3]]
mTable  = data.frame(ENSG = rownames(mTable), mTable)
AnnT = DESeq_Cohort1$AnnT
head(mTable)
                                                 ENSG     baseMean
ENSMUSG00000090025 Gm16088 ENSMUSG00000090025 Gm16088 2.013630e-01
ENSMUSG00000025900 Rp1         ENSMUSG00000025900 Rp1 7.937009e+00
ENSMUSG00000025902 Sox17     ENSMUSG00000025902 Sox17 2.084278e-02
ENSMUSG00000098104 Gm6085   ENSMUSG00000098104 Gm6085 2.618361e+00
ENSMUSG00000088000 Gm25493 ENSMUSG00000088000 Gm25493 5.372398e-02
ENSMUSG00000033845 Mrpl15   ENSMUSG00000033845 Mrpl15 1.402656e+03
                           log2FoldChange      lfcSE       stat
ENSMUSG00000090025 Gm16088   -0.028940390 0.09136916 -0.3167413
ENSMUSG00000025900 Rp1        0.106469263 0.14046002  0.7580040
ENSMUSG00000025902 Sox17      0.008723823 0.03637123  0.2398551
ENSMUSG00000098104 Gm6085     0.291415163 0.14016518  2.0790839
ENSMUSG00000088000 Gm25493    0.013212903 0.03637421  0.3632492
ENSMUSG00000033845 Mrpl15     0.056124421 0.02751339  2.0398950
                               pvalue      padj
ENSMUSG00000090025 Gm16088 0.75143989 0.8628992
ENSMUSG00000025900 Rp1     0.44844856 0.6361525
ENSMUSG00000025902 Sox17   0.81044263        NA
ENSMUSG00000098104 Gm6085  0.03760964 0.1171989
ENSMUSG00000088000 Gm25493 0.71641875        NA
ENSMUSG00000033845 Mrpl15  0.04136079 0.1256043
cat(nrow(filter(mTable, padj < 0.01, log2FoldChange > 0)), " up-regulated genes.\n", 
  nrow(filter(mTable, padj < 0.01, log2FoldChange < 0)), " down-regulated genes.\n", 
  sep = "")
1592 up-regulated genes.
1652 down-regulated genes.
hTable = mTable[rownames(AnnT)[AnnT[, "Homolog_Human"] != ""], ]
cat(nrow(filter(hTable, padj < 0.01, log2FoldChange > 0)), 
  " up-regulated genes with human homologs.\n", 
  nrow(filter(hTable, padj < 0.01, log2FoldChange < 0)), 
  " down-regulated genes with human homologs.\n", sep = "")
1201 up-regulated genes with human homologs.
1361 down-regulated genes with human homologs.

Pathway enrichment

source("/sibcb2/bioinformatics/Script/PathwayEnrichment.R")
load("DataSummary/pwDB.RData")
pwDB <- pwDB[c(2:5,9,14:18)]

UP = unique(AnnT[filter(hTable, log2FoldChange > 0, padj < 0.01)$ENSG, "Homolog_Human"])
DN = unique(AnnT[filter(hTable, log2FoldChange < 0, padj < 0.01)$ENSG, "Homolog_Human"])
PG = unique(AnnT[hTable$ENSG, "Homolog_Human"])

pwDBList <- list()
for(i in 1:length(pwDB))
  pwDBList[[i]] <- pwListFilter(pwDB[[i]], PG, minN = 15, maxN = 2000)
names(pwDBList) = names(pwDB)

write.table(as.matrix(pHyperListDB(UP, pwDBList, PG, FALSE)), 
  file = "out/rainbow_UP_PathwayEnrichment.txt",
  row.names = F, col.names = T, sep = "\t")
write.table(as.matrix(pHyperListDB(DN, pwDBList, PG, FALSE)), 
  file = "out/rainbow_DN_PathwayEnrichment.txt",
  row.names = F, col.names = T, sep = "\t")

cat(length(UP), " unique up-regulated genes with human homologs.\n",
  length(DN), " unique down-regulated genes with human homologs\n", sep = "")
1200 unique up-regulated genes with human homologs.
1361 unique down-regulated genes with human homologs
uTable = read.table(file = "out/rainbow_UP_PathwayEnrichment.txt", row.names = NULL, header = TRUE, sep = "\t")
uTable = arrange(uTable, Pvalue)
uPathways = c("HALLMARK_MTORC1_SIGNALING",
      "HALLMARK_OXIDATIVE_PHOSPHORYLATION",
      "HALLMARK_MYC_TARGETS_V1",
      "REACTOME_CHOLESTEROL_BIOSYNTHESIS",
      "REACTOME_REGULATION_OF_MITOTIC_CELL_CYCLE",
      "REACTOME_SIGNALING_BY_WNT",
      "REACTOME_METABOLISM_OF_MRNA",
      "REACTOME_ACTIVATION_OF_NF_KAPPAB_IN_B_CELLS")
uTable[uTable$PathwayName %in% uPathways, 1:7]
                                    PathwayName              DataBase
8                     HALLMARK_MTORC1_SIGNALING       MSigDB_h_all_v5
14           HALLMARK_OXIDATIVE_PHOSPHORYLATION       MSigDB_h_all_v5
17                      HALLMARK_MYC_TARGETS_V1       MSigDB_h_all_v5
24            REACTOME_CHOLESTEROL_BIOSYNTHESIS MSigDB_c2_cp_reactome
94    REACTOME_REGULATION_OF_MITOTIC_CELL_CYCLE MSigDB_c2_cp_reactome
97  REACTOME_ACTIVATION_OF_NF_KAPPAB_IN_B_CELLS MSigDB_c2_cp_reactome
111                   REACTOME_SIGNALING_BY_WNT MSigDB_c2_cp_reactome
134                 REACTOME_METABOLISM_OF_MRNA MSigDB_c2_cp_reactome
    SignatureSize PathwaySize Overlap       Pvalue       Qvalue
8            1200         196      61 5.957365e-22 4.557384e-19
14           1200         190      54 1.196260e-17 5.229363e-15
17           1200         193      52 5.675148e-16 2.043053e-13
24           1200          20      15 9.558724e-15 2.437475e-12
94           1200          76      24 4.603195e-10 2.996974e-08
97           1200          61      21 6.883535e-10 4.343014e-08
111          1200          61      20 4.458776e-09 2.458352e-07
134          1200         188      38 2.843030e-08 1.293631e-06
dTable = read.table(file = "out/rainbow_DN_PathwayEnrichment.txt", row.names = NULL, header = TRUE, sep = "\t")
dTable = arrange(dTable, -Pvalue)
dPathways = c("HALLMARK_HEME_METABOLISM",
        "HALLMARK_INTERFERON_GAMMA_RESPONSE",
        "HALLMARK_INTERFERON_ALPHA_RESPONSE",
        "HALLMARK_IL6_JAK_STAT3_SIGNALING",
        "HALLMARK_KRAS_SIGNALING_UP")
dTable[dTable$PathwayName %in% dPathways, 1:7]
                            PathwayName        DataBase SignatureSize
5792         HALLMARK_KRAS_SIGNALING_UP MSigDB_h_all_v5          1361
5852   HALLMARK_IL6_JAK_STAT3_SIGNALING MSigDB_h_all_v5          1361
5863 HALLMARK_INTERFERON_ALPHA_RESPONSE MSigDB_h_all_v5          1361
5971 HALLMARK_INTERFERON_GAMMA_RESPONSE MSigDB_h_all_v5          1361
6119           HALLMARK_HEME_METABOLISM MSigDB_h_all_v5          1361
     PathwaySize Overlap       Pvalue       Qvalue
5792         188      29 1.695827e-03 3.154548e-02
5852          82      16 9.805450e-04 2.230831e-02
5863          81      16 8.477694e-04 2.010988e-02
5971         179      31 1.558219e-04 6.315430e-03
6119         180      71 3.090972e-29 9.458375e-26
uTable = uTable %>% mutate(group = "UP", LogFDR = -log10(Qvalue))
dTable = dTable %>% mutate(group = "DN", LogFDR = log10(Qvalue))
pwTable = rbind(uTable, dTable)
colorG = rep("NC", nrow(pwTable))
colorG[pwTable$LogFDR > -log10(0.05)] = "upSig"
colorG[pwTable$LogFDR <  log10(0.05)] = "dnSig"
pwTable$colorG = colorG
pwTable = filter(pwTable, DataBase %in% c("MSigDB_h_all_v5", "MSigDB_c2_cp_reactome", "MSigDB_c2_cp_kegg", "MSigDB_c2_cp_biocarta"))

p <- ggplot(pwTable, aes(x = 1:nrow(pwTable), y = LogFDR, color = colorG))
p <- p + theme_bw() + labs(x = "Pathways", y = "-log10 FDR", title = "")
p <- p + theme(plot.title = element_text(hjust = 0.5)) + geom_point(size = 0.1, color = "grey")
p <- p + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p <- p + geom_point(size = 0.2)
p <- p + scale_colour_manual(values = c("upSig" = "#dd1c77", "dnSig" = "#3182bd", "NC" = "grey"))
print(p)

Overlapping CRISPR hits

load("DataSummary/Focused_Depleted.RData")

upGS = unique(sapply(strsplit(filter(mTable, log2FoldChange > 0, padj < 0.01)$ENSG, " "), function(z) z[2]))
dnGS = unique(sapply(strsplit(filter(mTable, log2FoldChange < 0, padj < 0.01)$ENSG, " "), function(z) z[2]))
pgGS = unique(sapply(strsplit(mTable$ENSG, " "), function(z) z[2]))
BG   = intersect(pgGS, PG)
UP   = intersect(upGS, BG)
DN   = intersect(dnGS, BG)

### MA cell line
res <- matrix(NA, length(maList), 6)
for(i in 1:length(maList)){

  Depleted = intersect(maList[[i]], BG)
  U1N = length(intersect(UP, Depleted))
  U0N = length(setdiff(UP, Depleted))
  D1N = length(intersect(DN, Depleted))
  D0N = length(setdiff(DN, Depleted))

  ABCD = matrix(c(U1N, U0N, D1N, D0N), nrow = 2, ncol = 2, byrow = TRUE)
  dimnames(ABCD) = list(c("UP", "DN"), c("Depleted", "NotDepleted"))

  res[i, 1:4] = c(U1N, U0N, D1N, D0N)
  res[i, 6] = fisher.test(ABCD)$p.value
  res[i, 5] = fisher.test(ABCD)$estimate
}
colnames(res) = c("UP_Depleted", "UP_NotDepleted", "DN_Depleted", "DN_NotDepleted", "OR", "Pvalue")
rownames(res) = names(maList)
knitr::kable(data.frame(res))
UP_Depleted UP_NotDepleted DN_Depleted DN_NotDepleted OR Pvalue
Post_vs_plasmid 102 23 39 22 2.4882466 0.0107318
Post_vs_Pre 92 33 29 32 3.0559849 0.0005989
Post_vs_vitro 15 110 13 48 0.5054908 0.1255157
Pre_vs_plasmid 100 25 31 30 3.8390473 0.0000710
vitro_vs_plasmid 98 27 31 30 3.4860765 0.0001872
vitro_vs_Pre 93 32 29 32 3.1847263 0.0004950
### HM cell line
res <- matrix(NA, length(hmList), 6)
for(i in 1:length(hmList)){

  Depleted = intersect(hmList[[i]], BG)
  U1N = length(intersect(UP, Depleted))
  U0N = length(setdiff(UP, Depleted))
  D1N = length(intersect(DN, Depleted))
  D0N = length(setdiff(DN, Depleted))

  ABCD = matrix(c(U1N, U0N, D1N, D0N), nrow = 2, ncol = 2, byrow = TRUE)
  dimnames(ABCD) = list(c("UP", "DN"), c("Depleted", "NotDepleted"))

  res[i, 1:4] = c(U1N, U0N, D1N, D0N)
  res[i, 6] = fisher.test(ABCD)$p.value
  res[i, 5] = fisher.test(ABCD)$estimate
}
colnames(res) = c("UP_Depleted", "UP_NotDepleted", "DN_Depleted", "DN_NotDepleted", "OR", "Pvalue")
rownames(res) = names(hmList)
knitr::kable(data.frame(res))
UP_Depleted UP_NotDepleted DN_Depleted DN_NotDepleted OR Pvalue
Post_vs_plasmid 102 23 36 25 3.0589979 0.0013226
Post_vs_Pre 99 26 38 23 2.2935265 0.0205806
Post_vs_vitro 12 113 8 53 0.7049096 0.4607979
Pre_vs_plasmid 95 30 25 36 4.5176767 0.0000042
vitro_vs_plasmid 95 30 37 24 2.0456464 0.0388821
vitro_vs_Pre 94 31 31 30 2.9158955 0.0014460

sessionInfo

R version 4.1.0 (2021-05-18)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server 7.6 (Maipo)

Matrix products: default
BLAS/LAPACK: /sibcb2/bioinformatics/software/Miniconda3/lib/libopenblasp-r0.3.15.so

locale:
[1] C

attached base packages:
[1] parallel  grid      stats     graphics  grDevices utils    
[7] datasets  methods   base     

other attached packages:
 [1] viper_1.26.0        survcomp_1.42.0     prodlim_2019.11.13 
 [4] gridExtra_2.3       survHD_0.99.1       survC1_1.0-3       
 [7] Hmisc_4.5-0         Formula_1.2-4       lattice_0.20-44    
[10] penalized_0.9-52    survival_3.2-11     Biobase_2.52.0     
[13] BiocGenerics_0.38.0 stringr_1.5.0       dplyr_1.1.2        
[16] scales_1.2.0        RColorBrewer_1.1-3  pheatmap_1.0.12    
[19] VennDiagram_1.6.20  futile.logger_1.4.3 ggrepel_0.9.1      
[22] ggplot2_3.4.2       rmarkdown_2.20     

loaded via a namespace (and not attached):
 [1] segmented_1.3-4      survivalROC_1.0.3    nlme_3.1-152        
 [4] tools_4.1.0          backports_1.2.1      bslib_0.2.5.1       
 [7] utf8_1.2.2           R6_2.5.1             KernSmooth_2.23-20  
[10] rpart_4.1-15         rmeta_3.0            DBI_1.1.3           
[13] mgcv_1.8-36          colorspace_2.0-3     nnet_7.3-16         
[16] withr_2.5.0          tidyselect_1.2.1     downlit_0.4.2       
[19] compiler_4.1.0       cli_3.6.1            formatR_1.11        
[22] htmlTable_2.2.1      labeling_0.4.2       bookdown_0.33       
[25] sass_0.4.0           checkmate_2.0.0      proxy_0.4-27        
[28] digest_0.6.29        mixtools_1.2.0       foreign_0.8-81      
[31] base64enc_0.1-3      jpeg_0.1-9           pkgconfig_2.0.3     
[34] htmltools_0.5.2      parallelly_1.32.0    fastmap_1.1.0       
[37] highr_0.9            htmlwidgets_1.5.4    rlang_1.1.1         
[40] rstudioapi_0.15.0    SuppDists_1.1-9.7    jquerylib_0.1.4     
[43] farver_2.1.1         generics_0.1.3       jsonlite_1.8.7      
[46] distill_1.5          magrittr_2.0.3       Matrix_1.3-4        
[49] Rcpp_1.0.9           munsell_0.5.0        fansi_1.0.3         
[52] lifecycle_1.0.3      stringi_1.7.8        yaml_2.3.5          
[55] MASS_7.3-54          listenv_0.8.0        splines_4.1.0       
[58] knitr_1.42           pillar_1.9.0         future.apply_1.9.0  
[61] codetools_0.2-18     futile.options_1.0.1 glue_1.6.2          
[64] evaluate_0.20        latticeExtra_0.6-29  lambda.r_1.2.4      
[67] data.table_1.14.2    png_0.1-7            vctrs_0.6.3         
[70] bootstrap_2019.6     gtable_0.3.0         kernlab_0.9-29      
[73] future_1.26.1        cachem_1.0.6         xfun_0.37           
[76] e1071_1.7-11         class_7.3-19         tibble_3.2.1        
[79] memoise_2.0.1        cluster_2.1.2        lava_1.6.10         
[82] globals_0.15.1